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Abstract: The phase reduction method is a dimension reduction method for weakly driven 
limit-cycle oscillators, which has played an important role in the theoretical analysis of synchro¬ 
nization phenomena. Recently, we proposed a generalization of the phase reduction method 
[W. Kurebayashi et al, Phys. Rev. Lett. Ill, 2013]. This generalized phase reduction method 
can robustly predict the dynamics of strongly driven oscillators, for which the conventional 
phase reduction method fails. In this generalized method, the external input to the oscillator 
should be properly decomposed into a slowly varying component and remaining weak fluctua¬ 
tions. In this paper, we propose a simple criterion for timescale decomposition of the external 
input, which gives accurate prediction of the phase dynamics and enables us to systematically 
apply the generalized phase reduction method to a general class of limit-cycle oscillators. The 
validity of the criterion is conhrmed by numerical simulations. 
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1. Introduction 

Synchronization of limit-cycle oscillators is a ubiquitous phenomenon that has been widely studied 
in many disciplines including physics, chemistry, biology, mechanical engineering, and electrical en¬ 
gineering [1-3]. The phase reduction method [1] has played a key role in the theoretical analysis 
of synchronization phenomena of weakly driven limit-cycle oscillators. This method enables us to 
reduce the dynamical equation of a high-dimensional limit-cycle oscillator to a one-dimensional phase 
equation, which facilitates theoretical analysis of various synchronization phenomena. Recently, en¬ 
gineering applications of the phase reduction method, e.g., dynamical analysis and optimal design of 
circuit and other oscillators [4-9] and optimal control of periodically spiking neurons [10-12], have 
been actively studied. 

However, the conventional phase reduction method has a drawback in practical applications, he., it 
works only when the external input to the oscillator can be assumed sufficiently weak. This limitation 
significantly narrows the applicability of the conventional method, because the weakness of the input 
cannot be assumed in many practical applications. In order to overcome this limitation, we recently 



proposed a generalized phase reduction method [13], which robustly works even for largely varying 
inputs under appropriate conditions. Using the generalized method, we can theoretically analyze the 
phase dynamics of strongly driven limit-cycle oscillators. 

When we use the generalized phase reduction method, we need to decompose the external input 
to the oscillator into a slowly varying low-frequency component and sufficiently weak fluctuations. 
Though this timescale decomposition can significantly affect the accuracy of the resulting phase 
equation, Ref. [13] did not provide an explicit criterion for decomposing the external input. In 
the present study, we propose a simple practical criterion to decompose the external input into low- 
frequency and high-frequency components, which yields a reasonable approximation of the oscillator 
dynamics. The proposed decomposition method will enable us to systematically apply the generalized 
phase reduction method to the analysis of various synchronization phenomena. 


2. Generalized phase reduction method 

We consider a limit-cycle oscillator driven by a general input I{t) — [di(t),... that 

smoothly depends on time t, described by 

^ = F(X(t),/(f)), (1) 

where X = [Xi ,..., X„]~'~ G K” is the state of the oscillator and F{X, I) ~ [Fi(X, /),..., Fn{X, /)]''' 
G M" is a vector field that represents the dynamics of the oscillator. We assume that there exists a 
finite interval A c K™ of the input value I such that the vector field F{X,I) has a stable periodic 
orbit Xo(t,/) with period T{I) and frequency uj{I) := 27r/T{I) when the input I is kept constant, 
and that this periodic orbit smoothly depends on / G A. 

Using the generalized phase reduction method [13], we can reduce the high-dimensional dynamics 
of the limit-cycle oscillator described by Eq. (1) to a one-dimensional phase equation. We define a 
generalized asymptotic phase Q{X,I) of the limit cycle XQ{t,I) that satisfies 


de{X,I) 


F{X,I) 


a;(/). 


( 2 ) 


for each constant I E A. We then decompose the input I{t) into a low-frequency component q{et) = 
[qi{€t), ..., G A and a high-frequency component ap{t) = a[pi{t),... ,Pm{t)]^ G K™ as 


I{t) = q{et) + o-p(t). 


(3) 


where the low-frequency component q(et) is assumed to vary slowly as compared to the amplitude 
relaxation time of the oscillator, and the high-frequency component crp{t) is assumed to be sufficiently 
weak. The small parameters e and a represent the slow timescale of the low-frequency component 
q{et) and the intensity of the high-frequency component ap{t), respectively. We introduce a phase 
variable 0{t) representing the state of the oscillator as 


eit) = eix{t),qiet)), 


(4) 


which depends on the slow low-frequency component q{et) of the input. We also define a conventional 
phase sensitivity function Z{6, q) G M" and two other sensitivity functions ^(0, q) G M™, q) G 
as follows: 


Z{e,q) 

i{0,q) 


dQ{X,q) 


dX 

dQ{X,q) 


dq 


X=Xo(e/c.(q),q) 


X=Xo{e/Lj{q),q) 


GiX,q)^Z{e,q)\^_ 


X=Xo{e/ui(q)..q) ’ 


(5) 

( 6 ) 
(7) 


where the (j, fc)-th element of G{X,q) G is given by G^^'^\X,q) := Then, the 

dynamics of the phase variable 6{t) is described by the following generalized phase equation [13]: 










dO{t) 

dt 


= U}{q{et)) + e^{e, q{et)) ■ q{et) + aC{0, q{et)) ■ p{t) 


O 


ea 


ea 


\{q{et)Y ’ \{q{et)) ’ X{q{et)) ’ \{q{et)Y 


( 8 ) 


Here, q{t) denotes dq{et)/d{et) and the function A(/) is the absolute value of the second largest 
Floquet exponent of the oscillator (1) driven by a constant input / G H, which characterizes the 
amplitude relaxation time of the oscillator. The generalized phase equation (8) is valid when 


A(g(et)) 


<C 1 and 


X{q{etW 


< 1 , 


(9) 


i.e., when the amplitude relaxation is sufficiently fast (see Ref. [13] for a detailed discussion). We 
hereafter assume that these conditions are satisfied. 

3. Simple criterion for timescale decomposition of external inpnts 

In the generalized phase reduction method, we need to decompose the input I{t) as in Eq. (3). How 
to decompose the input I{t) is an important problem, which can significantly affect the accuracy of 
the generalized phase equation (8). Our aim in this paper is to propose a simple criterion for choosing 
the threshold frequency 0^ that gives a reasonable decomposition of the input into low-frequency and 
high-frequency components for approximating the dynamics of the oscillator. In our derivation, the 
essential parameter for the timescale decomposition is the amplitude relaxation time of the oscillator; 
the statistical property of the external input I{t) (e.g., the power spectrum) is not important. 

We define the decomposition of the input I{t) by a linear filter /(r) as follows: 


/ -l-oo 

- T)f{T)dT, 

'OO 

ap{t) = lit)- qiet), 


( 10 ) 

( 11 ) 


where /(r) is assumed to be an ideal low-pass filter with the cutoff frequency fid, i-e., its amplitude 
response H(H) := | /(T)e“*^'^dT| of /(r) is given by 


A(fl) = 




i\n\ < Hd), 
0 (otherwise). 


( 12 ) 


As discussed in Appendix A, we can describe the dynamics of a limit-cycle oscillator by the phase 
and amplitude variables, where the amplitude variable represents the deviation of the oscillator state 
from the periodic orbit. In particular, when the oscillator state AC(t) is two-dimensional, it can be 
fully described by a phase variable 0(t) defined in Eq. (4) and an amplitude variable r(t) defined as 


r{t) = R{Xit),qiet)), 
where the function A(X, J) of X G R" and I G R™ satisfies 

dR{X,I) 


dx 


FiX,I) = -XiI)RiX,I). 


(13) 


(14) 


As shown in Appendix A, we can derive the following dynamical equation for the amplitude variable 
r(f): 


dr{t) 

dt 


= -A(q(et))r -h aCr{0, r, q{et)) ■ pit) + e^riO, r, qiet)) ■ qiet) -f 0(o-^), 


(15) 


where = GiX,q) 


T dRjX.q) I 


iriO,r,q\ 


_ dR{X,q) I 


, and X(0, r, q) is a 


dX IX=X(S,r-.q)’ I 5 dq IX=X(0,r,q)’ 

state point in R" satisfying 0(X,g) = 6 and RiX,q) = r. This equation shows that the amplitude 
r(t) fluctuates around r = 0 due to the external input. 














The approximation error of the generalized phase equation (8) is 0{r) (see Appendix A). Thus, we 
can minimize the approximation error by minimizing the deviation r{t) from the periodic orbit. As 
shown in Appendix B, we can approximate 0, 1) by Cr{0, 0, 1) as follows: 


Uo,o,i) 


1 

a(7) 


Cr{o,o,i) + o 



Moreover, from Eq. (15), the order of r can be evaluated as follows: 


r = O 


_ 


By plugging Eqs. (16) and (17) into Eq. (15), we can obtain 


dr{t) 

dt 


-X{q{et))r + Cr{0, 0, q{et)) ■ I{t) + O 


\X{q{et)) 


ea \ 

A(g(et))’ A(g(et))y ' 


(16) 


(17) 


(18) 


where I{t) is a transformed external input, whose j-th element is given by 

ij{t) := apj{t) + (19) 

for j = 1,..., m. We define the variance of Ij{t) as 

Vj{nd) = lim - [ [Ij{t)]'^dt, (20) 

T^oo r Jq 

which is finite because I{t) is smooth and bounded. For smaller Vj{fld), the fluctuation of r{t) 
becomes smaller and the generalized phase equation will give more precise prediction of the phase 
dynamics. 

To derive a simple criterion for determining the threshold frequency 12^, we assume that the decay 
rate X{q{et)) of r{t) does not vary too violently and thus its typical values can be characterized by 

Ac := X{q{et)) ~ ^ “ J q{^t)d^ , (21) 


where q{et) is the long-time average of the slowly varying part of the input, q{et). Though this is a 
rather rough characterization of the decay rate of r{t), it enables us to derive a simple criterion for 
the threshold frequency. Replacing X{q{et)) in Eq. (19) with Ac, Eq. (20) can be estimated as follows: 

r^d q 2 I'oa 

Rj(fld)w2/ ^Pj{n)dn + 2 Pj{n)dn, (22) 

Jo X^ 

where Pj{Xl) is the power spectrum of Ij{t). The optimal threshold frequency Xld = 12^ that minimizes 
this approximate variance Vj{Xld) can be determined as 


Xl*d = Xr 


(23) 


because this 12^ satisfies 


r^d / q 2 \ 

v,[nd)-v,{n*d) = 2 

r^d / q 2 \ 

v,{nd)-v,{n*d) = 2 


Pj{nd)dn > 0 , 
Pj{nd)dn > 0, 


for ild<^*d, 
for Qd > XId- 


(24) 

(25) 


Thus, under the above approximating assumptions, the optimal timescale for the decomposition of 
the input that minimizes the variance Vj{Xld) of the input coincides with the characteristic amplitude 
relaxation time of the oscillator. 














Fig. 1. Time series of the input I{t) to the oscillator for j = 5, and the 
low-frequency and high-frequency components g(£t) and ap(t) decomposed by 
low-pass filters with the threshold frequencies fid = 0.5, 2, 5 and 10. 


We propose Eq. (23) as a simple criterion for choosing the value of the threshold frequency fid- It 
gives the optimal fid for predicting the oscillator dynamics when A(g(et)) is strictly constant, and 
is expected to provide a reasonable prediction even if A(g(et)) varies slowly. The criterion (23) is 
valid for general external inputs, including periodic signals with delta-peaked power spectra, because 
we did not introduce any assumptions on the statistical property of the input in the derivation 
of the criterion (23). Also, though we assumed that the state of the oscillator is two-dimensional 
for simplicity, the above result can be generalized to higher-dimensional cases by regarding A(g) as 
the absolute value of the second largest Floquet exponent among the n Floquet exponents of the 
oscillator, because the deviation from the periodic orbit is dominated by the slowest amplitude mode 
characterized by the second largest Floquet exponent. 


4. Numerical simulations_ 

In order to confirm the validity of our criterion (23), we performed numerical simulations for evaluating 
the effect of the threshold frequency fid on the approximation accuracy of the generalized phase 
equation (8). In the simulation, we use a modified Stuart-Landau oscillator [1] dehned as 

^ = F,(x,y,I(t)) := e2^W(x - y - I) - [(x - I(t)f + y\x - /(f)), (26) 

= F2{x,y,I{t)) := /) - [(x - I{t)f +y'^]y, (27) 

where X = [x, is a state variable, and /(f) is the input to the oscillator that is decomposed 

into a low-frequency component q{et) and a high-frequency component crp{t). For this oscillator, the 

absolute value of the second largest Floquet exponent is A{q) = 2e^'^. 

Defining the phase 0{t) for this oscillator as 













































Fig. 2. Mean square errors of the generalized phase equation ( 8 ) versus the 
threshold frequency fid for 7 = 5, 7 and 10. The arrow represents the threshold 
frequency fid = 2 given by the criterion (23). 


6{t) = 0{x{t), y{t), q{et)) = tan’ 


vit) 


x{t) - q{et) ’ 

we can reduce Eqs. (26) and (27) to the following generalized phase equation: 
dOit) 


dt 


= sin 6 » ■ q{et) + cosO] -pit). 


We generated the input to the oscillator I(t) by a Fourier series given by 


500 


I{t) = sm{aet + Pe), 


£=1 


where Pi is an i.i.d. random number drawn from a uniform distribution in [0,27r], 


Ai = JlO-^ 


1 + ’ 
q,^ = 0.1(£-0.5), 


(28) 

(29) 

(30) 

(31) 

(32) 


and 7 is a parameter representing the characteristic scale of I{t). In this case, the long-time average 
of I{t) becomes zero. Thus, the criterion (23) gives a threshold frequency 


Ac = A(</)|^^^ = 2. (33) 

We computed the time series of the state variable X{t) — [x{t),y{t)]^ through direct numerical 
simulations of Eqs. (26) and (27), and evaluated the approximation accuracy of the generalized phase 
equation ( 8 ) by the mean square error MSE(f7d) defined as 


MSE(nd) = 


7' 


^exact 


(t) - e, 


approx 


(7 


dt, 


(34) 


where 0 approx(i) is a predicted value of 0{t) obtained by plugging 9{t) = Q{x{t),y{t),q{et)) into the 
generalized phase equation ( 8 ), and 0exact(i) is the exact value of 9{t) given by 


77 = 


de{x{t),y{t),q{et)) 

dx 


Fi{x{t),y{t),I{t)) -h 


dQ{x{t),y{t),q{et)) 

dy 


F 2 {x{t),y(t),I(t)) 


de{x{t),y{t),q{et)) dq{et) 
dq 


dt 


(35) 















which can be directly calculated from the time series of X{t), q{et) and I{t). Note that the exact 
value of 0{t) (Eq. (35)) also depends on the value of fid, because the definition of the phase variable 
0{t) (Eq. (34)) itself depends on 

Figure 1 shows the time series of I{t), q{et) and up{t) for = 0.5, 2,5 and 10 and 7 = 5. Figure 2 
shows the mean square error MSE(n(;). We see that MSE(r2,^) has a minimum for each value of the 
parameter 7 , which is in reasonable agreement with the criterion (23), i.e., Ac = 2. Though the 
criterion (23) proposed in this paper is based on a simplifying assumptions, these results indicate that 
our criterion is able to give a reasonable threshold frequency 0 ^. 


5. Conclusion 

In this paper, we proposed a simple criterion of the threshold frequency for timescale decomposition 
of the external input for generalized phase reduction of limit-cycle oscillators. Under the assumptions 
that the amplitude relaxation is sufficiently fast and the timescale of the amplitude relaxation can be 
characterized by the second largest Floquet exponent of the oscillator when it is driven by a constant 
long-time average of the input, we derived a criterion for choosing the threshold frequency, which 
is simple and physically reasonable. We confirmed the validity of our criterion by direct numerical 
simulations. The criterion proposed in this paper is simple and easy to use, and thus it will be helpful 
in the engineering applications of the generalized phase reduction method, e.g., optimal design of 
circuit oscillators [7,8] and optimal control of periodically spiking neurons [10,11]. 
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Appendix A: Variable transformation 

In the following appendices, we briefly review the derivation of the generalized phase equation (8), 
including the transformation of the state variable X (t) to the phase variable 0{t) and the amplitude 
variable r{t), and the relation between the sensitivity functions, Eq. (16). The results shown here are 
the same as those given in the Supplementary Information of our previous paper, Ref. [13]. 

We consider a limit-cycle oscillator whose dynamics depends on a time-varying input I{t): 

X{t) = F{Xit),Iit)). (36) 


The state variable X{t) is assumed to be two-dimensional here, but the result can be extended to 
higher-dimensional cases. As argued in the Supplementary Information of Ref. [14] by Goldobin et 
ah, we can define a phase 9 = Q{X,I) and an amplitude r = R{X,I) of the oscillator satisfying 


de{X,I) 

dR{X,I) 

ax 


F(X,7) 

F{X,I) 


w(J), 

-X{I)R{X,I), 


(37) 

(38) 


where A(J) is the absolute value of the second Floquet exponent of the oscillator for constant I. Thus, 


9{t) = uj{I), r{t) = —X{I)r 


(39) 


when the input I is constant and in the given range A c R™. 

When the parameter I{t) varies with time, we decompose I{t) into a slowly varying component 
q{et) and remaining weak fluctuations (7p{t) as I{t) = q{et) + crp{t), and define the phase 9{t) and 
the amplitude r{t) of the oscillator as 


9{t) = e{X{t),q{et)), 
r{t) = R{X{t),q{et)). 

The dynamical equations for d{t) and r{t) are then given by 


9 = 


r = 


deix,i) 


dX 

dR{X,I) 


dX{t) dQ{X,I) 


dX 


(X,q(et)) 


(X.q(et)) 


dt 


dl 


dX{t) dR{X,I) 


dt 


dl 


(X,q{et)) 


iX,q{et)) 


dq{et) 

dt 

dq{et) 

dt 


Plugging I{t) = q{et) + (jp{t) into Eq. (36) and expanding it in a, we obtain 

X = F(X, q{et)) + aG{X, q{et))p{t) + 


(40) 

(41) 


(42) 

(43) 

(44) 


where G(X, q) is the matrix defined in the main article as G^^’^\X, q) := (j, k = 1,2 here). 

Substituting Eqs. (37), (38), and (44) into Eqs. (42) and (43), we obtain 


9 = uj{q{et)) + ( 


de{X,I) 


dX 


G{X,q{et))p{t) +e 


dQ{X,I) 


(X.q{et)) 


dl 


q{et)+ 0{a‘^), (45) 


{X,q(et)) 


f = -X{q{et))r + i 


dR{X,I) 


dX 


■ G{X,q{et))p{t) -h e 


dR{X,I) 


(X.q(£t)) 


dl 


(X,q(£t)) 


• qiet) + 0 ( 0 - 2 ), 


(46) 


where qiet) = dqiet)/diet). By defining Q{9,r,I) G M™, (^r{9,r,I) G R™, ^g{9,r,I) G R™ and 
^ri9,r,I) G R™, respectively, as 


























X=X{e,r,I) 


(47) 


CeiO,r, I) 


CriO,r,I) 

ie{0,r,I) 

UO,r,I) 


G{X, 

G(X, 


IV 

IV 


de{x,i) 

M 

dR{X,I) 

ax 


x=x(e,r,i) 


de{x,i) 

al 

dR{X,I) 

a7 


X=X{e,r,I) 


X=X(e,r,I) 


(48) 

(49) 

(50) 


where X{0,r,I) e represents an oscillator state with 0 = Q{X,I), r = R{X,I), and parameter 
I, Eqs. (45) and (46) can be written as 


0 = iLi{q{et)) + aCe{0,r,q{et)) ■ p{t) + e^g{0,r,q{et)) ■ q{et) + 0{cr^), (51) 

r = -\{q{et))r + (jQr{0,r,q{et))-p{t) + eir{0,r,q{et))-q{et) + O{a‘^). (52) 


Here, and ^@(0,0,/) with r = 0 correspond to the sensitivity functions C(^)-f) ^{0,1) 

defined in the main article. The other two functions Cr{0,f,I) and ^r{0','f’,I) represent sensitivities 
of the amplitude variable to the small fluctuations and to the slowly varying component of the input, 
respectively. 

When X{q{et)) is sufficiently large, the amplitude r{t) takes tiny values around 0, and the phase 
0{t) approximately obeys Eq. (51) with r = 0, i.e., the generalized phase equation (8), which yields 
an approximation error of 0{r). See Ref. [13] for a detailed discussion on the validity of the approxi¬ 
mation. 


Appendix B: Derivation of Eq. (16) 

In this Appendix, we derive Eq. (16). As explained in the main article, we assume the existence of 
a stable limit-cycle orbit Xo(t,/) with frequency u){I) satisfying dX^it, I)/dt = F{XQ{t, I), I) that 
smoothly depends on I for each constant IgA. We denote the phase of the oscillator as 0{t) = u>{I)t, 
and represent the oscillator state as Xo{t,I) = Xo{0{t),I) by using the phase 0{t) instead of t. 

We differentiate both sides of Eq. (38) with respect to I and plug in X = Xq{0,I). From the 
left-hand side, we obtain 


d 



\— 

fdR{X,I)\-] 

Fix, I) 

dl 

[ dX J 

x=Xoiej) 

[dX 

[ dl J\ 



C.(0,O,/), (53) 


x=Xoie,i) 


where we used the definition of G{X, I) and Eq. (48) with r = 0. The first term on the right-hand 
side can be further calculated as 


d /dR{X,I) 

'M V 


dl 


Fix, I) 


= UJ(I) 


\ ^ 

fdRiX,I)\] 

dX 

[ dl )\ 


x=Xo{e.i) 

x=Xo{e,i) 


\ ^ 

^ai?(x,/)Ai 

dX 

^ 91 )\ 


dXo{0it),I) 


x=Xo{e,i) 

dXoi0,I) d f dRjX,!) 

d0 ^^^’d0 V dl 


dt 


X=Xo{B,I) 


= u}(r 


dU0,0,I) 

d0 


(54) 


where we used the chain rule for the derivative d/d0 and Eq. (50). Similarly, by differentiating the 
right-hand side of Eq. (38) with respect to I, we can derive 


— [-A(I)R(X,J) 


X=Xa(B.I) 


d\{I) 

dl 


i?(X,J) + A(J) 


dRiX,I) 


X=Xo{B,I) 


= -A(J)^,(0,O,J), 
(55) 


where we used R{X, = 0. Thus, we obtain a linear first-order ordinary differential 

equation for 4r-(^, 0, J), 










































+ Cr{0, 0, I) = -A(/)C,(0, 0, I), 


which can be solved as 
1 


Uo,o,i) = - 


e-Xr 0 


_ — [ (gi Q ^ _J_ f 

By expanding the integrand, the order of ^r{S,0X) can be estimated as 

UO, 0 ’ ^) = 0 , 1)ds + o 

Cr{0,0,I)+O 


HI) 


HI) 


which gives Eq. (16). 










